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We propose a method to study dynamical response of a quantum system by evolving it with an 
imaginary-time dependent Hamiltonian. The leading non-adiabatic response of the system driven to 
a quantum-critical point is universal and characterized by the same exponents in real and imaginary 
time. For a linear quench protocol, the fidelity susceptibility and the geometric tensor naturally 
emerge in the response functions. Beyond linear response, we extend the finite-size scaling theory 
of quantum phase transitions to non-equilibrium setups. This allows, e.g., for studies of quantum 
phase transitions in systems of fixed finite size by monitoring expectation values as a function 
of the quench velocity. Non-equilibrium imaginary-time dynamics is also amenable to quantum 
Monte Carlo (QMC) simulations, with a scheme that we introduce here and apply to quenches of 
the transverse-field Ising model to quantum-critical points in one and two dimensions. The QMC 
method is generic and can be applied to a wide range of models and non-equilibrium setups. 



I. INTRODUCTION 

The dynamics of thermally isolated quantum systems 
beyond linear response has become a focus of experimen- 
tal and theoretical research in thermalization^i^ univer- 
sal quantum critical dynamics?^ quantum annealing;^ 
and many other areasi^ It has been realized^ that de- 
viations from adiabaticity in gapless systems and near 
quantum-critical points, in particular, can be charac- 
terized by scaling behavior of the fidelity susceptibility 
and its adiabatic generalizations. These susceptibilities 
are related to non-equal time correlations of the cor- 
responding quench operator—!^ evaluated either at the 
beginning or the end of the dynamical process. One 
can, thus, extract valuable information on the dynam- 
ical properties of quantum systems by analyzing their 
non-adiabatic response. Such response can be directly 
measured experimentalljiiSiii or studied numerically. At 
the moment, numerical studies of real-time dynamics of 
interacting systems are limited to small systems, mostly 
in one dimension, however 

We here show that quantum dynamics can also be sim- 
ulated by evolving the system in imaginary time. In par- 
ticular, we demonstrate that the leading non-adiabatic 
response of a system with its Hamiltonian changing in 
imaginary time is very similar to that of the real-time dy- 
namics. This allows us to use powerful quantum Monte 
Carlo (QMC) techniques to investigate the dynamical re- 
sponse. Another advance presented here is the extension 
of the standard linear response theory, both in real and 
imaginary time, to non-linear driving protocols where the 
velocity or acceleration of the quench replaces the ampli- 
tude. In particular, we show that the linear response of 
physical observables in the case of linear quenches is char- 
acterized by the components of the geometric tensorj^ii^ 
This allows one to experimentally measure them, or simu- 
late them using QMC, and study their singularities near 
quantum critical points. Previously, with QMC simu- 
lations it was only known how to compute the diago- 
nal elements of the geometric tensor, i.e., the fidelity 
susceptibilities 



We show that the non-perturbative response of generic 
observables can be described by extending the standard 
finite-size scaling theory of quantum phase transitions to 
non-equilibrium protocols (e.g., by simultaneous scaling 
in the system size and the quench velocity, or by only 
changing the velocity at fixed system size), with expo- 
nents that we derive here. In this work we focus on 
imaginary time dynamics, but all universal results also 
apply to real-time protocols. 

We discuss the underlying time-evolution formalism 
and results of adiabatic perturbation theory in Sec. HIl 
followed by results from linear response theory of phys- 
ical observables to the quench velocity and the emer- 
gence from it of the geometric tensor in Sec. IIIII Then, 
in Sec. IIVI we formulate the scaling theory of non- 
perturbative response of interacting systems to slow per- 
turbations near quantum-critical points, extending the 
scaling theory of phase transitions to non-equilibrium 
setups. In Sec. |V] we apply the theory to the particu- 
lar example of the one-dimensional transverse-field Ising 
model, where the scaling forms derived can be compared 
with exact results. In Sec. I VII we present the QMC 
method and numerical results obtained with it for the 
two-dimensional transverse-field Ising model. Wc con- 
clude in Sec. IVIII with a brief summary and discussion. 
More details of the adiabatic perturbation theory are 
given in Appendix]^ and properties of the dynamic sus- 
ceptibilities derived are further discussed in Appendix IB] 



II. TIME EVOLUTION 

Wc consider the imaginary-time evolution described by 
a Hamiltonian 'H(A) which implicitly depends on time 
through the tuning parameter A(t). We assume that the 
evolution starts at some time tq < and ends at r = 0, 
A(0) being the point of interest. To simplify the notation 
we set A(0) = 0. The imaginary-time propagation of the 
wave function in this setup is governed by the Shrodingcr 
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equation in imaginary time: 



dMr) = -7{(A(t))^(t). 



(1) 



The formal solution at time r is given by the evolution 
operator, iP{t) = U'iP{tq), with U given by the time- 
ordered exponential: 



U — Tr cxp 



dT'n{x{T')) 



(2) 



Before going into details of the dynamics, let us make 
some remarks: (i) In the adiabatic limit, A — > 0, the 
system rapidly falls into its instantaneous ground state, 
after a transient time, and then follows this state, (ii) At 
finite A the system is constantly excited from the ground 
state by the evolving Hamiltonian and relaxes back due 
to imaginary-time propagation. The proximity to the in- 
stantaneous ground state is controlled by A(r) near the 
final point t ^ 0. If the velocity vanishes at the final 
point, A(0) = 0, then the degree of nonadiabaticity is 
controlled by the acceleration A(0), etc. (iii) Imaginary 
time evolution is amenable to QMC simulations, giving 
access to universal aspects of quantum dynamics in a 
wide range of systems. We will outline such a gener- 
alization of standard equilibrium QMC in Sec. I VII and 
apply it to the transverse-field Ising model. We first dis- 
cuss the analytical framework needed for analyzing both 
QMC and experimental results. 

The easiest way to analyze the general properties of the 
solution of Eq. ([!]) is to go to the adiabatic (co-moving) 
basis. This procedure is similar to that in real time, 
though containing very important subtleties. The de- 
tails of the analysis are given in Appendix Here we 
present only the final result of the first order of adiabatic 
perturbation theory, which contains all relevant scaling 
information. Denoting by a„ (0) the expansion coefficient 
of the wave-function ip{0) in the eigenstates of the final 
Hamiltonian we have 



a„(0) 



dT 



{n\dr'H\Q) 
A«o(t) 



exp 



dr' A„o(r') 



(3) 



where /^^^{t) = f„(T) — £q{t) is the instantaneous en- 
ergy of the n-th level relative to the ground state and 
(ri|9^H|0) is the transition matrix element between the 
instantaneous eigenstates. 

To make further progress in analyzing the transition 
amplitudes ([3|), let us consider the very slow asymptotic 
limit A — 0. To be specific, we assume that near t ~ 
the tuning parameter has the form A(r) k, v\T^\/r\ 
(see also Ref. Q). The parameter w, which controls the 
adiabaticity, plays the role of the quench amplitude (r = 
0), velocity (r = 1), acceleration (r = 2) etc. It is easy to 
check that in the asymptotic limit u ^ 0, Eq. ([3]) gives 



{n\dxn\Q) 



(4) 



where all matrix elements and energy levels are evaluated 
at r = 0. From this perturbative result we can evaluate 
the leading non-adiabatic response of various observables 
and define the corresponding susceptibilities. 



III. LINEAR RESPONSE AND GEOMETRIC 
TENSOR 

Let us represent the observables of interest as gener- 
alized forces, i.e., derivatives of H with respect to the 
couplings fj,] = —dfj^H. By using this representation 
we do not lose any generality. For example, a spin-spin 
correlation function • of some lattice model can be 
represented as a response with respect to an infinitesimal 
coupling connecting these spins. Then we find 



where C = (i/'(0)|Mp|V'(0)) and 

1 ^ {0\^xn\n){n\^^H\0)+^i^\ 



(r+l) 



2{£n - fo)'-+i 



(5) 



(6) 



These susceptibilities can also be expressed through the 
imaginaryiS and real time connected correlation functions. 



(r+l) 
v - 



1 

21? 



(7) 

where OxHt is the imaginary-time Heisenberg representa- 
tion of the operator d\}-i evaluated at t (and in real time 
one substitutes r ^ it -I- 0"*", as discussed in Appendix 
IB|) . Thus, by changing the exponent r of the quench pro- 
tocol one can probe different moments of the real and 
imaginary time correlation functions. Let us point out 
that the factor L~'^ in Eqs. (jG]) and (O is inserted for 
convenience for extensive observables, which appear as a 
response to global perturbations. For intensive observ- 
ables this factor is not needed. 

The situation is slightly different for diagonal observ- 
ables like the energy. 



Q = {n)- £o 



the log-fidelity," 
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F = -ln(|(^(0)|0)n, 



(8) 



(9) 



or the entropy entropy which in the lowest order of per- 
turbation theory are described by quadratic rather then 
linear responses^ 



Q 



2 (2r+l) 

« Xaa ' 

r ^ 2 (2r+2) 



(10) 



The response coefficients in Eq. ([5]) have a very interest- 
ing geometric interpretation for linear quenches (r = 1). 
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Then the susceptibihty xl^x reduces to the symmetrized 
/iA-component of the geometric tensorf^iH which defines 
the Ricmannian metric in the manifold of the ground 
states of the Hamiltonian A).— The diagonal com- 

(2) 

ponents of the geometric tensor X\\ define the fidelity 
susceptibilities!^ We emphasize that the metric ten- 
sor, which was originally thought to have no physical 
significance^ emerges here as a response of physical ob- 
servablcs to the quench velocity. Thus, Eq. ^ opens a 
practical (numerical or experimental) way of analyzing 
the geometry of the ground state wave function in the 
parameter space and studying its universality, nature of 
its singularities, and its topology. 



IV. SCALING THEORY 

In gapped systems all non-equal time correlation func- 
tions decay exponentially with time, implying that the 
susceptibilities x^™^ converge for all m in the thermody- 
namic limit. For gapless systems the situation is more 
complicated and the susceptibilities can diverge. To un- 
derstand the nature of this divergences we will employ 
scaling analysis. 

If the quench operator dxH is marginal or relevant, 
then its scaling dimension is. 



Aa s dim[9A^] = 2 — dim[A]. 



(11) 



For a marginal perturbation maintaining gaplessness in 
the vicinity of A = and not affecting the dynamic ex- 
ponent 2, we have dim [A] = and Ax = z. Such cases 
include superfluids and Fermi liquids (with A the inter- 
action coupling). If dxH is relevant, e.g., when driving 
the system to a gapped phase at A ^ 0, then by def- 
inition dim[A] = l/j^,where v is the correlation length 
exponentiii Then = z — 1/i^, and from Eq. ([5]) we 
obtain 

r,'-;^!''^ ^ dim[xi'r = A, + d - 1/^ - zr. (12) 
For jj, = X this expression reduces to a known result 



If 77 



(r+1) 



tcm size, 



< the susceptibility diverges with the sys- 



(r+l) 



(13) 



and the perturbative result ([S]) breaks down in the ther- 
modynamic limit. To find the correct asymptotics of the 
observables in this case, we introduce the scaling dimen- 
sion of the velocity; dim[u] = dim[A + zr = l/u -\- zr. 
Following arguments similar to Ref. H instead of Eqs. ^ 
and pop we then find: 



(14) 



Here C is some non-universal constant. However, unlike 
in Eq. ([5]), this constant docs not have to be the ground 



state expectation value. Whether the ground state ex- 
pectation value is included or not in C determines the 
small velocity asymptotics of the scaling functions /^a (x) 
and ffix{x). In general these asymptotics can be deter- 
mined from physical arguments. For x ^ I we should 
recover linear or quadratic response for diagonal and off- 
diagonal observables, respectively, plus possibly universal 
ground state contribution if it is not included in C. The 
large argument asymptotic of the scaling function can 
be obtained from other considerations. For example, for 
extensive operators the scaling functions should saturate 
at large x so that is extensive. The properties of the 
susceptibilities ([6]) are further discussed in Appendix [Bj 
Instead of the length scale equal to the system size L 
in Eq. ([H)) there can be another relevant length scale, 
e.g., the distance X12 = |xi — X2[ between two points 
xi and X2 if we are interested in correlation functions. 
Thus, in a translationally invariant system one expects 
that the non-equilibrium connected correlation function 
in a large system should scale as: 

(Af^(xi)M4x2))e « -^f {vxi;^'^-) . (15) 

2^12 

Likewise we can generalize Eq. ([T?|) to quenches which at 
the final time end up in the vicinity of the QCP, i.e., at 
A/ A,: 



vV 



' „l/(zi/r+l) I • ^^^> 



This scaling relation can be used for independently locat- 
ing the quantum critical point by, e.g., sweeping across 
the phase transition in a sufficiently big system with dif- 
ferent velocities (such that vL'^^^l^ ^ 1). Then there 
will be a crossing point in M^/?;('^+^'')'^/(^+^^) plotted 
versus the coupling A/ in curves corresponding to differ- 
ent velocities. The expression (fT6|) . which is applicable 
to real-time protocols as well, also suggests a convenient 
way for determining the location of the critical point ex- 
perimentally, by changing the velocity for a fixed sys- 
tem size. The usual finite-size scaling procedure requires 
changing the system size, which is not always feasible. 
On the other hand, changing the quench velocity would 
normally be quite straightforward in experiments on, e.g., 
cold atoms. 

The scaling relations generalize the standard 

finite-size scaling theory of quantum phase transitions 
(corresponding to r = 0) to non-equilibrium protocols 
and constitute our main analytical result. Eq. (|14p is 
valid both for real and imaginary time, and also for the 
expectation value (il/^) (in which case the small argu- 
ment asymptotics is dictated by the requirement that 
(M^) reduces to its equilibrium value in the adiabatic 
limit), as well for any other observable O with scaling 
dimension A©. For example, for the particular cases of 
the energy and the fidelity, A^; = 2 and A^;^ = 0, respec- 
tively, and Eq. (jl4[) reduces to known results;^ which were 
recently verified numerically for a particular ID model 
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Observable 




Q 


F 


vL^ < 1 




-128-*' ^ 




vL^ > 1 


0.26V^L 


0.0265 vL 


O.O2760yL 



TABLE I: Scaling of the excess interaction energy Ez with 
respect to the final ground state, the excess total energy Q, 
and the log-fidelity F with the quench rate and the system 
size for the transverse-field Ising chain. 

V. ID TRANVERSE-FIELD ISING MODEL 

To illustrate the above general scaling results we 
consider a linear quench in the one-dimensional (ID) 
transverse-field Ising model with Hamiltonian 



i (ij) 



where ax and Uz are the Pauli matrices and (ij) are near- 
est neighbours sites. The dimensionless coupling con- 
stant J drives the system through a critical point at 
Jc = 1, with critical exponents z = v = 1. We consider 
the following quench protocol: J(t) = 1 + A = 1 + vt, 
starting in the ground state at tq = —1/v. Using 
the Jordan- Wigner transformation, the model can be 
mapped to free fermions. The analysis of the imagi- 
nary time dynamics is straightforward and available in 
the literature for similar real-time setupsi ^^i^° We there- 
fore only quote our results in Table ID 

It is evident that the general scaling prediction 
indeed applies to this example. To illustrate further the 
finite-size scaling behavior predicted by Eq. ([Ti|) . in Fig.[T] 
(left panel) we plot the shift of the interaction energy 
with respect to the final ground state. 



= -J 



(ij) 



(ij) 



|0) 



(18) 



versus vL"^. The data for different system sizes collapse, 
showing that one can use the proposed imaginary-time 
adiabatic approach to extract critical properties of a 
quantum phase transition. To test our predictions fur- 
ther, we analyze the square of the longitudinal magneti- 
zation (the order parameter); 



(19) 



which has scaling dimension = l/4»ii This together 
with Eq. dHI) imply that 




(20) 



The large and small argument asymptotics of the scaling 
function f^j are dictated by the equilibrium asymptotics 
in the diabatic limit, fzj{x) ~ const at x ^ 1, and by the 
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FIG. 1: (Color online) Excess interaction energy E^ (left) and 
squared magnetization (right) of Ising chains graphed accord- 
ing to our scaling predictions. The lines show the expected 
asymptotic forms; the slope is 1 for vL^ <^ 1 and 1/2 for 



(17) ^'^ > 1 on the left and -3/8 for vL > 1 on the right. 



requirement that ^ 1/L at vL^ 3> 1 when quenching 
from the disordered phase. If we quench from the ordered 
phase, J(to) > 1, then fzj{x) ~ x^^^ at a; 1 (so 
that ml const). The finite-size scaling predictions and 
asymptotics are in excellent agreement with numerical 
data (Fig. [1] left panel) obtained using QMC simulations 
with the algorithm discussed next. 



VI. QUANTUM MONTE CARLO METHOD 

A major advantage of the imaginary-time approach is 
that generalized QMC methods can be applied to evolve a 
state with the operator ^ . Here we use an approach sim- 
ilar to the stochastic series expansion (SSE) method, dis- 
cussed in the context of the transverse-field Ising model 
in Ref. [2l|. The method is generally applicable to all mod- 
els for which standard equilibrium QMC simulations can 
be used, i.e., those for which there is no sign problem. 
Below we first briefly review standard finite-temperature 
and ground-state QMC approaches. We then outline 
the general idea of the non-equilibrium QMC (NEQMC) 
method in imaginary time and apply it to the one- and 
two-dimensional transverse-field Ising models. 



A. Standard QMC methods 

Standard QMC algorithms can be classified into finite- 
temperature methods, where the goal is to compute a 
quantum-mechanical thermal average of the form 



(A) 



Tr{e-^-^} 



(21) 



and ground-state projector methods, where some oper- 
ator P(/3) is applied to a "trial state" |^'o), such that 
'5/3) = P(/3)|^'o) approaches the ground state when 
/? — > oo. Normally one is interested in expectation values. 



(A) 



1 



i^Ml'^p), Z={^p\^p), (22) 
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which approache the corresponding true ground state ex- 
pectation values, (A) — > (0|A|0}, when /3 — cx). For 
the projector, one can use the imaginary-time evolution 
operator P(/3) = = e'^'^, with a fixed (time- 

independent) Hamiltonian, or one can use a high power of 
the Hamiltonian, P(/3) = where /? oc m/N gives the 
same rate of convergence (which is governed by the gap 
between the ground state and the first excited state in 
the symmetry sector of the trial state) for the two choices 
for a given system volume N. This follows from a Taylor 
expansion of the time evolution operator, which for large 
P is dominated by powers of the order n = /3|i?o|, where 
Eq is the ground state energy (and Eq oc N) . 

There are several ways to deal with the exponential. 
In the context of spins and bosons, the most frequently 
used methods are based on (i) the Suzuki- Trotter- 
decomposition, which leads to world-line methods,— 
(ii) the continuous-time version of world-lines (e.g., the 
worm algorithro^i) and (iii) the Taylor expansion lead- 
ing to the SSE method2ii^i2£ (see Ref. [13 for a recent 
review of these approaches) . The latter two methods are 
not affected by any approximations (beyond statistical 
sampling errors), while (i) has a discretization error. 




FIG. 2: (Color online) Sampled imaginary-time sequences 
shown versus the propagation index p after two successive 
Monte Carlo sweeps in a linear-quench simulation of an 8 x 8 
2D transverse-field Ising model. At the initial time tq — the 
Hamiltonian contains only the transverse field {J = 0, h — 1). 
The J-term is increased linearly with time to the critical point 
(J « 3.05, /i = 1) at the final time t — 2. There is a total 
of 2m = 1280 operators, with p = 1, . . . , 640 in a term of the 
projection of the ket state in ((23| and p — 641, . . . , 1280 in 
a corresponding bra term. The set of time points shown in 
red was obtained from those shown in black by the method 
of updating several overlapping segments of a large number 
of times (here approximately 100) as discussed in the text. 
The right panel shows the behavior close to the center of the 
string (the final time) in greater detail. 



B. Non-equilibrium QMC algorithm 

The NEQMC algorithm is similar to a ground-state 
projection, but instead of e~^'^ for a fixed Hamiltonian 
one uses the evolution operator ^ with a time dependent 
Hamiltonian. As in equilibrium QMC, one can treat the 
exponential operator in several different ways. Here we 
employ the series expansion. 

Evolving from tq to r, Eq. ^ is expanded in a power- 
series and applied to an initial state |^'(0)): 

[~niT,,)]-.-[~n{nmm. (23) 



Writing —Ti in terms of individual site and bond opera- 
tors, here denoted Hi, i = 1, . 



(24) 



the operator product is written as a sum over all strings 
of these operators. Truncating at some maximum power 
n = m (adapted to cause no detectable truncation error, 
as in the SSE method^l) and introducing a trivial unit 
operator Hq = 1, we can write Eq. as 



l*(r)) 



■^-^ (m — n)! 

iJ,,„(T,„)---i?,,(T2)iJ,,(Tl)|*(0)), 









■ dT2 dn X 




J To J To 



(25) 



where ip G {0, . . . , A'op}, X]_f/ ^he sum over all se- 
quences ii, . . . and n is the number of indices ip ^ 



in a given sequence. More generally, beyond the 
transverse-field Ising model, i would refer to a lattice unit 
as well as a diagonal or off-diagonal part of the operator 
on this unit. The operators Hi then have have the prop- 
erty that Hi\a) = hi{a)\a'), where \a') is a basis state, 
i.e., in the basis chosen to expand the states, there is no 
branching of the series of states obtained in the sequence 
of states resulting from the operators acting one-by-one 
in Eq. (^5]). 

As always in QMC simulations, we are in practicve 
restricted to systems for which the expansion is positive- 
definite, which is the same class for which sign problems 
can be avoided in equilibrium simulations. While the 
sign problem is a limitation of the QMC approach in 
general, the class of accessible models is still large and 
includes highly non-trivial and important systems. With 
the series expansion used in the NEQMC method here, 
avoiding the sign problem places constraints on the ma- 
trix elements ft.,; (a) — the product of all matrix elements 
corresponding to a term in has to be positive. 

Expectation values 



{A)r 



{^{r)\A\^>{r)) 



(26) 



are computed by sampling the normalization (^'(t) |\I'(t)) 
written with (f25|) . For the transverse-field Ising model, 
which we will apply the method to below, the method is 
very similar to the one developed in Ref. [2l| in the context 
of SSE QMC, the main difference being the change in 
the time boundaries; from periodic at finite temperature 
to those dictated by the initial state |vE'(0)} of the time 
evolution. Changes in the operator sequence are made 
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with the times fixed. The times are updated separately. 

The operator sampling in the case of the transverse- 
field Ising model is particularly simple when the starting 
state is the equal superposition, 



N 



l*(0)) 



(27) 



which we use below, but other states can be used as well 
(in particular, for Heisenberg and other spin-isotropic 
systems, amplitude-product states in the valence-bond 
basis^^ are very convenient, and a generalization of the 
loop updates used in the ground-state projector method 
of Ref. 12^ can be used). 

Since efficient operator and state cluster-updates have 
been described in detail in the literature for various 
models in standard QMC simulations ;^^'^^ only the time 
update (which is a generalization of a scheme previ- 
ously discussed for equilibrium QMC in the interaction 
representation^H) will be briefly outlined here. A whole 
segment of times, r^, . . . , r^+n, can be simultaneously up- 
dated by generating n -I- 1 numbers within the range 
(Ti_i , Ti+„+i), then order these times according to a stan- 
dard scheme scaling as log(ri))^ and inserting the ordered 
set in place of the old segment of times. The Metropolis 
acceptance probability is easily obtained from (j25p . at 
a cost scaling as n. The number n can be adjusted to 
give an acceptance probability close to 1/2. Fig. [2] shows 
an example of a time sequence and how it changes af- 
ter a sweep of updates of partially overlapping segments 
covering the whole sequence of times. 



C. Results 

Using the NEQMC method we first confirmed that the 
exact results for the Ising chain are reproduced. Com- 
plete agreement was found to within very small statistical 
errors. The results in the right panel of Fig. [1] are from 
the NEQMC simulations. 

We next considered the same model on the 2D square 
lattice, i.e., the generalization of the ID Hamiltonian 
(|T7l) . The critical coupling in this case is Jc = 0.32841 
(based on exact diagonalization of a series of small lat- 
tices, which show behavior agreeing very well with pre- 
dictions from low-energy field theory) ^ In the left panel 
of Fig. |3] we show the scaling of the excess Ising energy 
Ez, i.e., the 2D generalization of p^ . for i x L lattices 
with L up to 64, using the known^ exponent ly — 0.6298 
(obtained for the classical 3D Ising model, which should 
be in the same universality as the 2D quantum model 
studied here, which has dynamic exponent when z = 1). 
We have divided by the leading powers of L and v 
predicted above and, hence, we should obtain a constant 
behavior for large x. This is not quite seen yet for these 
systems sizes, but the eventual convergence seems plau- 
sible. For smaller x the data collapse very well and the 
asymptotic x — > behavior is reproduced. In the right 
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FIG. 3: Same as Fig.[T]for the square-lattice model (quenching 
to its critical point Jc ~ 0.32841). QMC data are scaled 
according to the theoretical results, with the lines illustrating 
the predicted asymptotic slopes with exponents z — 1, v = 
0.6298, and r; = 0.0364 (note that = l + ri)M 



panel we show that also the squared magnetization scales 
according to our predictions, over five decades of the scal- 
ing argument t;L*^''+^-'/^. 



VII. SUMMARY AND DISCUSSION 



We have shown that detailed information on static and 
dynamic properties of a system can be obtained by prop- 
agating it in imaginary time. There are many similarities 
with real-time dynamics. In particular, we showed that 
one can use imaginary time to obtain universal exponents 
characterizing quantum critical points and to measure 
the fidelity susceptibility and components of the geomet- 
ric tensor as response of physical observables to a linear 
quench. We obtained finite-size scaling expressions char- 
acterizing the response of various observables with the 
quench rate. In this way we extended the scaling theory 
of quantum phase transitions to non-equilibrium proto- 
cols. A clear advantage of the imaginary-time approach 
is that one can use powerful QMC simulations and cir- 
cumvent complications related to real-time simulations. 
We have presented such a generic non-equilibrium QMC 
scheme and illustrated this approach for the transverse- 
field Ising model. Exact results (in one dimension) and 
QMC results (in one and two dimensions) show excel- 
lent agreement with the scaling predictions. The QMC 
method will be useful for studies of a wide range of non- 
trivial models on large lattices. 

The ideas presented here apply also to quantum an- 
nealing, i.e., protocols where, in order to analyze the 
ground state of a complicated classical or quantum prob- 
lem, one introduces an auxiliary coupling which makes 
the Hamiltonian simple and then slowly decreases this 
coupling to zero. This will allow one to address quantum 
annealing problems using QMC simulationsi^ 
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where A„o(t) = £n{T)—£o{T). The matrix element above 
for non-degenerate states can also be expressed as: 

{n\dr\0) = -(n|9,H(r)|0)/A„o(r). (A7) 

Appendix B: Adiabatic susceptibilities and 
non-equal time correlation functions 



Appendix A: Adiabatic perturbation theory 

Let us discuss the leading non-adiabatic correction to 
the imaginary-time Schrodinger equation ([T]): 

dMr) = ~n{x{T))^. (Ai) 

The natural way to address this question is to use adia- 
batic perturbation theory (APT), similar to that devel- 
oped in Refs. [20ll35| in real time. We write the wave func- 
tion in the instantaneous eigenbasis {|7i(A))} of H{X): 

V'(r) = ^a„(r)|n(A(r))). (A2) 

n 

Substituting this expansion into Eq. ([T]) wc find 

^ +^amiT)(n\dr\m) = -£„(A)a„(T), (A3) 

where £'„(A) are the eigenenergics of H(A) corresponding 
to the states \n). Making the transformation 

rO 



a„(T) = a„(T)exp 



£n{T')dT' 



(A4) 



we can rewrite Eq. ([T]) as an integral equation [and note 
that Q!„(0) = a„(0)]: 

q;„(t) = a„(0) -I- ^ / dr' {n\dr'\m)ajn{T') 



X exp 



dr" [£„(r")~f™(r")] 



(.A5) 



In principle one should supply this equation with ini- 
tial conditions at r = tq but, as we argued earlier, it is 
not necessary if |ro| is sufficiently large, since the sensi- 
tivity to the initial condition will be exponentially sup- 
pressed. Instead we impose the asymptotic condition 
a„(T — ^ — oo) — >■ Sno, implying that far in the past the 
system is effectively in the ground state. 

Eq. (|A5p is convenient for analysis with the APT. In 
particular, if the rate of change is very small, A(t) 0, 
then to leading order in A the system remains in the 
ground state; am{T) « Smo (except during the initial 
transient, which is unimportant at large |to|). In the 
next higher order the transition amplitudes to the states 
n ^ are given by: 



a„(0) 



dr {?i\dr\0) exp 



dr' A„o(r') 



(A6) 



In this appendix we discuss the properties of the adia- 
batic susceptibilities [Eq. © of the main text] : 



(r+l) 



{0\dxH\n){n\d^H\Q) + ^ X 
2(£„ - £oy+' 



(Bl) 



For linear quenches these quantities reduce to the sym- 
metrized components of the geometric tensor— up to a 
normalization factor. The representation of these suscep- 
tibilities through imaginary time correlation functions is 
a straightforward generalization of the result contained 
in Ref.i (see also Ref.i): 



(r+l) 



where 



1 

21? 



dxUr = c'^dx'Hc-^'^. 



(B2) 



(B3) 



Performing the Wick's rotation t ^ it -\- where e is 
an infinitesimal positive number, we extend this result to 
real time: 



(r+l)^ 

X^x 



2L<^ 



dt*-{0\d^ntdxHo + dxntd,,no\0),, 



where 



dx-Ht 



dxHc 



-itn 



(B4) 



(B5) 



stands for the real-time Heiscnberg operator. Thus we 
see that the adiabatic susceptibilities of order r+l probe 
the r-th moment of the symmetric retarded correlation 
function of the operators dxTi. and d^J-L. Introducing the 
Fourier transform of this correlation function: 



dte'^'{0\dt,'Htdx-Ho + dxHtd^Ho\0)c, 

(B6) 



we see that the susceptibility X^'a''^'' '^^^ expressed 
through derivatives of the imaginary part of functions 



which define the structure factors: 



(r+l) 

V 



1 



ImG:,^ '[lj) 



fiX 



(B7) 



Finally let us mention the representation of these sus- 
ceptibilities through the real part of non-equal time cor- 
relation functions. This can be achieved cither by apply- 
ing Kramers-Kronig relations to the equation above or 
directly from the definition: 
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(r+l) 



"5^0 



{Q\dx'H\n){n\d^n\Q) + 11 ^ \ 



5{£n - So 



did 



(B8) 



where 

-oo 

(B9) 



1 T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 
900 (2006). 

2 M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 
(2008). 

^ A. Polkovnikov, Phys. Rev. B 72, 161201(R) (2005). 

* W. H. Zurek, U. Dorner, and P. ZoUer, Phys. Rev. Lett. 
95, 105701 (2005). 

^ A. Das and B. K. Chakrabarti (Editors), Quantum An- 
nealing and Related Optimization Methods, Lecture Note 
in Physics, vol. 679 (Springer, Heidelberg, 2005). 
J. Dziarmaga, Adv. in Phys. 59, 1063 (2010). 
A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat- 
tore, Rev. Mod. Phys. 83, 863 (2011). 

^ C. De Grandi, V. Gritsev, and A. Polkovnikov, Phys. Rev. 
B 81, 012303 (2010). 

^ L. C. Venuti and P. Zanardi, Pliys. Rev. Lett. 99, 095701 
(2007). 

1° C.-L. Hung, X. Zhang, N. Gemelke, and C. Chin, Phys. 

Rev. Lett. 104, 160403 (2010). 
" D. Chen, M. White, C. Borries, and B. DeMarco, 

arXiv:1103.4662 (2011). 

M. Kolodrubetz, D. Pekker, B. K. Clark, and K. Sengupta, 
arXiv:1106.4031 (2011). 

J. P. Provost and G. Vallee, Comm. Math. Phys. 76, 289 
(1980). 

A. F. Albuquerque, F. Alet, C. Sire, and S. Capponi, Phys. 
Rev. B 81, 064418 (2010). 

M. M. Rams and B. Damski, Phys. Rev. Lett. 106, 055701 
(2011). 

S.-J. Gu and H.-Q. Lin, Europhys. Lett. 87, 10003 (2009). 
S. Sachdev, Quantum Phase Transitions (Cambridge Uni- 



18 



versity Press, 1999). 
D. Schwandt, F. Alet, and S. Capponi, Phys. Rev. Lett. 



103, 170501 (2009). 

J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005). 

C. De Grandi and A. Polkovnikov, Led. Notes in Phys. 

802, 75 (2010). 

A. W. Sandvik, Phys. Rev. E 68, 056701 (2003). 

M. Suzuki, S. Miyashita and A. Kuroda, Prog. Theor. 

Phys. 58, 1377 (1977). 

J. E. Hirsch, R. L Sugar, D. J. Scalapino, and R. Blanken- 
becler, Phys. Rev. B 26, 5033 (1982). 

N. V. Prokof'ev, B. V. Svistunov, and I. S. Tupitsyn, Sov. 
Phys JETP 87, 310 (1998) [arXivxond- mat/9703200]. 
A. W. Sandvik and J. Kurkijarvi, Phys. Rev. B 43, 5950 
(1991). 

A. W. Sandvik.J. Phys. A 25, 3667 (1992). 

A. W. Sandvik, AIP Conf. Proc. 1297, 135 (2010) 

(arXiv:1101.3281). 

S. Liang, B. Doucot, and P. W. Anderson, Phys. Rev. Lett. 
61, 365 (1988). 

A. W. Sandvik and H. G. Evertz, Phys. Rev. B 82, 024407 
(2010). 

^° A. W. Sandvik, R. R. P. Singh, and D. K. Campbell Phys. 

Rev. B 56, 14510 (1997). 
^'^ W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. 

T. Vetterling, Numerical Recipes: The Art of Scientific 

Computing (Cambridge University Press, 2007). 

C. J. Hamer, J. Phys. A: Math. Gen. 33, 6683 (2000). 

M. Hasenbusch, K. Pinn, and S. Vinti, Pliys. Rev. B 59, 

11471 (1999). 

C.-W. Liu, A. Polkovnikov, A. W. Sandvik, work in 
progress. 

G. Rigolin, G. Ortiz, and V. H. Ponce, Phys. Rev. A 78, 
052508 (2008). 



